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ABSTRACT 

This  report  describes  a  new,  flexible  computer  code  in  the  FORTRAN  computer 
language  to  make  ray  calculations  for  ocean  acoustic  tomography.  The  Numerical 
Recipes  software  package  provided  the  basis  for  much  of  this  computer  code.  The  ray 
equations  are  reviewed,  and  ray  equations  that  include  the  effects  of  ocean  current  are 
derived.  Methods  are  derived  for  rapidly  integrating  those  equations  to  obtain  time  front 
and  eigenray  information  for  long-range,  deep-water  acoustic  transmissions.  These  meth¬ 
ods  include  a  look-up  table  for  sound  speed,  sound  speed  gradient,  second  derivative  of 
sound  speed,  and  range-dependent  information.  Cubic  spline  methods  are  used  to  inter¬ 
polate  sound  speed  with  depth  and  to  obtain  the  derivatives  of  sound  speed.  The  choice 
of  the  step  size  increments  used  to  integrate  the  equations  is  a  critical  aspect  of  the  inte¬ 
gration,  affecting  both  the  accuracy  of  the  prediction  and  the  speed  of  computation.  A 
predetermined,  user-specified  step  size  appears  to  allow  more  efficient  calculations  than 
"adaptive  step"  methods.  "Adaptive  step"  methods  adjust  the  step  size  automatically  to 
maintain  a  given  accuracy  in  the  integration  of  the  ray  equations,  while  user-specified 
step  sizes  allow  one  to  use  prior  knowledge  of  the  integration  problem  to  achieve  the 
desired  accuracy  with  much  less  computational  overhead.  Several  integration  methods 
were  explored,  but  the  classical  4th  order  Runge-Kutta  method  appears  to  be  the  most 
efficient  and  best  method  for  this  integration  problem.  Appendices  describe  detailed 
aspects  of  the  computer  code,  as  well  as  the  methods  used  for  deriving  eigenray  informa¬ 
tion  and  for  parallelizing  the  ray  calculations.  The  computer  code  is  designed  to  be 
unstable  so  that  the  user  can  easily  modify  it  to  his  or  her  own  porpoises. 
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MOTIVATION 

Although  numerous  ray  tracing  codes  are  available,  none  satisfy  all  of  the  present 
requirements  of  long-range  ocean  acoustic  tomography.  Ocean  acoustic  tomography  is 
described  by  Munk,  Worcester,  and  Wunsch  (1995).  What  is  required  for  long-range 
tomography  is  a  fast,  accurate,  and  flexible  code.  The  requirement  of  flexibility  necessi¬ 
tates  coding  that  can  be  easily  modified  by  the  user;  thus  FORTRAN  is  the  preferred 
computer  language.  Flexibility  allows  the  user  to  easily  implement  his  or  her  own  algo¬ 
rithms,  such  as  a  better  search  for  eigenrays  or  an  alternative  (perhaps  faster  or  more 
accurate)  integration  routine.  Issues  of  numerical  accuracy  in  ray  predictions  are  dis¬ 
cussed  in  the  section  on  integrating  the  ray  equations.  Finally,  the  requirement  of  speed 
necessitates  coding  that  can  sometimes  become  opaque.  This  report  therefore  also 
describes  the  methods  employed  to  achieve  computer  code  that  is  highly  efficient. 

The  present  and  immediate  goal  of  the  authors  of  this  report  is  to  achieve  fast, 
accurate  wavefront  and  eigenray  travel  time  predictions  at  basin  scale  ranges  (3—5  Mm) 
in  the  North  Pacific  Ocean  as  part  of  the  North  Pacific  Acoustic  Laboratory  and  Acoustic 
Thermometry  of  Ocean  Climate  projects  (Dushaw  1999;  Dushaw  et  al.  1999;  Worcester 
et  al.  1998,  Colosi  et  al.  1999).  These  projects  require  eigenray  predictions  for  a  half 
dozen  or  so  time  series  of  sound  speed  sections  derived  from  acoustic  data  obtained  at 
3-5  Mm  range.  For  a  single  sound  speed  section,  calculation  of  the  time  series  of  ray 
travel  times  takes  about  30  hours  on  a  200-MHz  Pentium  Pro  computer,  so  the  need  for 
the  fastest  possible  code  is  evident. 

Many  of  the  ideas  used  in  implementing  the  code  to  be  discussed  here  originated  in 
the  RAY  code  and  associated  technical  report  by  Bowlin  et  al.  (1992).  Alas,  this  code  is 
in  C.  The  new  computer  code  described  here  uses  the  FORTRAN  cubic  spline  and  inte¬ 
gration  routines  from  Numerical  Recipes  (Press  et  al.  1992;  Numerical  Recipes  here¬ 
inafter);  thus  use  of  this  code  requires  a  license  to  use  the  Numerical  Recipes  software 
($40). 

CREDO:  The  code  described  here  is  meant  to  be  easily  modified  by  the  user,  and 
so  it  will  never  be  a  stable  ray  code.  The  code  is  meant  to  be  transparent  and  fast. 
Metaphorically,  if  the  Bowlin  RAY  code  is  a  Fiat  2000  with  a  Fiat  engine,  the  code 
described  here  is  a  '67  Chevy  Impala  with  JATO  (Jet  Assisted  Take-Off)  propulsion.* 


♦Urban  Legend  No.  37 
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RAY  EQUATIONS,  OR  EQUATIONS  OF  MOTION 

According  to  Bowlin  et  al.  (1992):  The  equations  of  motion  for  a  ray  traveling 
through  the  ocean  can  be  cast  in  Cartesian  coordinates  as  follows: 


drc  „  dzc 

—  tan  8  -  — 

(la) 

c  c 

I 

II 

■SI 

(lb) 

dr 

dt  seed 

(lc) 

dr  c 

where  8  is  the  angle  of  the  ray  with  respect  to  the  horizontal  r  axis,  and  z  is  the  vertical 
coordinate  [positive  upward]. 

Bowlin  et  al.  (1992)  continue:  For  long  range  ocean  acoustics,  the  curvature  of  the 
Earth’s  surface  makes  non-Cartesian  coordinates  more  suitable  for  ray  tracing.  Let  new 
z’  axes  lie  along  radii  passing  through  the  center  of  the  Earth  with  z  —  0  at  sea  level  and 
z'  =  Re  at  the  Earth ’s  center,  where  Re  is  the  radius  of  the  Earth,  and  let  the  new  r  be  the 

range  measured  along  a  circular  arc  at  sea  level .  The  new  equations  of  motion  which 

include  the  geometrical  effects  due  to  a  spherical  Earth  are 
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These  are  the  equations  that  RAY  integrates....  [These  latter  equations  have  z  positive 
downward.]  However,  the  RAY  code  actually  integrates  not  0,  but  cos  6  and  sin  6,  thus 
avoiding  the  calculation  of  the  transcendental  functions. 

An  alternate  way  to  integrate  these  equations  is  to  use  the  equations  in  Cartesian 
coordinates  (the  equations  above  with  Re— »oo  and  fe  — >  1)  but  apply  the  well-known 
Earth  flattening  transformation  (Aid  and  Richards  1980).  This  transformation  is  applied 
once  to  the  initial  sound  speeds  and  associated  depths,  and  the  subsequent  integration  in 
the  Cartesian  coordinates  is  then  mathematically  equivalent  to  the  above  equations.  Phys¬ 
ically,  this  transformation  is  a  stretching  of  depth  and  sound  speed  equivalent  to  the  cur¬ 
vature  of  the  Earth’s  surface.  If  e  =  z/Re  (z  positive  downward)  then  the  Earth  flattening 
transformation  is  z  =  z  *  (1  +  e/2  +  e  *  e/3)  and  c  =  c  *  (1  +  e  +  e  *  e).  Since  some  com¬ 
putation  can  be  saved  if  this  transformation  is  applied  once  before  integrating  the  differ¬ 
ential  equations,  this  transformation  together  with  the  Cartesian  differential  equations  is 
preferred.  Ray  predictions  using  the  two  methods  agree  to  within  1  ms  at  3-Mm  range, 
and  the  integration  of  the  Cartesian  equations  appears  to  be  roughly  10%  faster. 

If  cs  =  cos  6,  sn  =  sin  6  and  we  apply  the  flat  Earth  transformation  to  z  and  c,  the 
equations  above  can  be  re-written  as 

dcs  fd7'C  9,'C  sn 

_  -  sn\ - 

dr  \  c  c  cs 


dsn  fdz'C  dr’C  sn 
dr'  c  c  cs 


(4b) 


dz  _  sn 
dr'  cs 


(4c) 


dt  _  1 

dr'  c  cs 


(4d) 


No  trigonometric  functions  need  to  be  calculated  while  integrating  these  equations.  Since 
integration  of  these  equations  is  most  sensitive  to  the  angle  integration,  the  redundant 
equations  for  angle  (the  equations  for  cs  and  sn)  are  actually  helpful  for  the  stability  of 
the  integration.  The  above  equations  can  be  cast  in  only  three  equations,  where  the  first 
equation  is  for  d(tan6)/dr,  but  this  formulation  involves  a  term  sqrt(l+tan  0*tan  6)  in  the 
integration  for  travel  time  (see  Appendix  B).  This  square  root  is  computationally  expen¬ 
sive  (for  one  test  case  it  increased  computation  time  by  about  15%),  perhaps  more  expen¬ 
sive  than  integrating  the  additional  differential  equation.  Similarly,  it  is  important  to  code 
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these  equations  using  the  smallest  number  of  divisions;  a  division  is  about  as  computa¬ 
tionally  expensive  as  taking  a  square  root. 

For  the  open-ocean  environment  of  the  North  Pacific,  the  term  involving  dc/dr  in 
the  above  equations  appears  to  be  5-6  orders  of  magnitude  less  than  the  other  terms.  Ray 
predictions  using  Levitus  sound  speeds  with  and  without  the  dc/dr  term  were  insignifi¬ 
cantly  different.  As  Bowlin  et  al.  (1992)  comment,  if...  an  environment  where  dropping 
the  dc/dr  term  [in  Eq.  4]  produces  a  significant  change  in  the  wavefront  then  it  is  very 
likely  that  the  linear  range  dependent  model  is  inadequate  for  that  environment  even  if 
the  dc/dr  term  is  included.  This  term  is  therefore  commented  out  in  the  code  for  effi¬ 
ciency  purposes,  but  the  user  concerned  with  significant  range  dependence  should  con¬ 
sider  checking  if  this  term  is  important  or  not. 


LOOK-UP  TABLES  AND  SOUND  SPEED  INTERPOLATION 

It  is  clear  that  a  fast  code  will  necessarily  rely  on  look-up  tables  for  sound  speed, 
sound  speed  gradient,  and  other  parameters.  Depth  is  the  obvious  index  to  these  look-up 
tables,  since  with  this  index  the  sound  speed  values  can  be  quickly  obtained  at  arbitrary 
depths  during  the  integration.  As  pointed  out  by  Bowlin  et  al.  (1992),  the  rapid,  accurate 
calculation  of  sound  speed,  sound  speed  gradient,  etc.,  at  arbitrary  depth  and  range  is  a 
critical  aspect  of  the  ray  prediction. 

For  the  present  code,  we  adopt  cubic  splines  for  interpolation  of  sound  speed 
{Numerical  Recipes )  and  for  the  calculation  of  sound  speed  gradients.  Splines  have  been 
discussed  in  the  context  of  calculating  acoustic  rays  by  Moler  and  Solomon  (1970). 
Cubic  splines  are  mathematically  the  most  natural  and  smoothest  interpolation;  their  goal 
in  life  is  to  minimize  the  second  derivative  of  the  interpolation  (Parker  1994).  The 
smoothest  possible  interpolation  is  essential  because  of  the  sensitivity  of  the  ray  predic¬ 
tions  to  sound  speed  gradient.  In  addition,  cubic  splines  allow  the  sound  speed  gradient 
to  be  calculated  quickly  and  accurately. 

We  note  again  that  it  is  essential  that  smoothly  varying  sound  speeds  be  used  for 
ray  predictions.  It  is  therefore  important  for  the  user  of  any  ray  trace  code  to  apply  any 
desired  smoothing  or  interpolation — horizontally  or  vertically — prior  to  the  ray  predic¬ 
tion.  The  code  does  not  implement  any  smoothing  of  the  user’s  sound  speed  profiles. 
Such  smoothing  necessarily  requires  allowing  some  misfit  to  the  sound  speed  values, 
which  can  lead  to  biases  in  the  sound  speed  profiles.  Since  the  code  described  by  Bowlin 
et  al.  (1992)  applies  some  smoothing  to  the  sound  speed  profiles,  Bowlin  et  al.  discuss 
this  bias  problem  at  length;  it  is  perhaps  a  problem  best  left  to  the  user  of  the  code. 

The  cubic  spline  interpolation  method  is  described  by  Numerical  Recipes  and 
Parker  (1994),  so  only  the  details  appropriate  to  the  present  application  are  discussed 
here.  Typically,  sound  speed  values  for  ray  tracing  are  available  at  a  small  number  of 
depths,  30-100  say.  The  cubic  spline  method  is  used  to  interpolate  those  values,  as  well 
as  values  for  the  first  and  second  derivatives,  to  3-m  increments  throughout  the  water 
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column,  and  these  values  are  used  to  fill  a  table  of  sound  speed  values.  (To  properly  han¬ 
dle  surface  reflected  rays,  the  value  for  sound  speed  at  the  surface  is  repeated  at  3  m 
above  the  ocean  surface.)  For  arbitrary  depth  z,  the  sound  speed  can  be  calculated  by 


c 


=  c  i 


(5) 


where  Cj  =  c(zj),  dz  =  z-  int(z),  and  j  =  int(z)  gives  the  index  in  the  sound  speed  look¬ 
up  table.  While  this  is  an  approximation  for  sound  speed,  it  is  sufficiently  accurate  for 
ray  tracing. 

Sound  speed  gradient  in  the  ray  equations  must  be  calculated  carefully,  since  even 
small  errors  in  this  quantity  result  in  ray  predictions  with  poor  quality.  With  cubic 
splines,  the  sound  speed  gradient  at  arbitrary  depth  z  between  the  y'th  and  (j  +  l)th  depths 
is  ( Numerical  Recipes ) 


dc  Cj+i  -  Cj 
dz  z  j+ 1  Zj 


(Zy+J  Zj)  C J  + 


(.Zj+i  Zj)Cj+ 1 


(6) 


where 


A  =  Zj+1  B  =  —  . =  1  -  A  (7) 

Z  j+ 1  Zj  Zj+ ]  Zj 

However,  the  difference  between  adjacent  depths  in  the  table  is  a  constant  Az  =  3  m,  and 
many  of  the  terms  in  this  equation  are  constant  and  so  may  be  calculated  prior  to  the  inte¬ 
gration.  When  a  1-m  depth  increment  is  used,  the  equations  simplify  even  further  at  the 
expense  of  a  very  large  look-up  table.  Thus,  the  equation  for  the  first  derivative  becomes 
(positive  z  is  downward) 


dc 

dz 


A2c " 
~iKz 


-  const. j 


(8) 


where  one  half  the  second  derivative  and  const. j  are  calculated  ahead  of  time  in  the  look¬ 
up  table. 


const,  j 


Cj+ 1  Cj  ^z{c  j  Cj+ j) 

Az  6 


(9) 


TM  3-98  5 


UNIVERSITY  OF  WASHINGTON  •  APPLIED  PHYSICS  LABORATORY 


The  vertical  gradient  of  sound  speed  can  be  calculated  rapidly  and  accurately  in  this  way. 
Note  that  while  the  Taylor  expansion  approximation  for  c(z )  in  Eq.  5  is  adequate,  the 
Taylor  expansion  approximation  for  dddz  is  not. 

Range  dependence  is  implemented  by  assuming  a  constant  horizontal  sound  speed 
gradient  and  a  constant  horizontal  gradient  of  the  vertical  gradient  of  sound  speed.  Thus, 


c(r,  z)  =  c{r j,  z)  +  (r  -  rj) 


c(jj+ 1,  z)-c(rj,  z) 
rM  ~  rJ 


(10) 


for  an  integration  step  at  range  r  between  the  y'th  and  (j  +  l)th  sound  speed  profiles.  The 
equation  for  the  horizontal  dependence  of  sound  speed  vertical  gradient  is  similar.  All 
indications  (Bowlin  et  al.  1992)  are  that  these  approximations  for  range  dependence  are 
adequate.  It  is  essential  to  implement  at  least  a  constant  sound  speed  gradient  between 
sound  speed  profiles,  however,  because  jumping  from  one  profile  to  another  in  a  discon¬ 
tinuous  fashion  produces  a  multitude  of  false  caustics  in  the  time  front. 

The  sound  speed  look-up  table  containing  the  information  required  to  calculate 
sound  speed  and  sound  speed  gradient,  including  their  range  dependence,  at  arbitrary 
depth  and  range  may  be  constructed  using  a  three-dimensional  matrix  of  depth,  range, 
and  six  variables  at  those  depths  and  ranges.  We  use  1 835  depths  at  3-m  increments  and 
the  ranges  given  by  the  initial  set  of  sound  speed  profiles.  The  six  variables  at  those 
depths  and  ranges  are  sound  speed,  sound  speed  gradient,  sound  speed  second  derivative 
with  depth,  horizontal  derivative  of  sound  speed,  horizontal  derivative  of  sound  speed 
gradient  (assuming  a  constant  gradient  between  profiles),  and  the  constant  term  described 
above  for  calculating  sound  speed  gradient  at  arbitrary  depth.  With  these  variables,  sound 
speed  and  sound  speed  gradient  can  be  calculated  at  arbitrary  depth  and  range  very  effi¬ 
ciently.  This  double  precision  table,  if  filled,  uses  about  13  Mbytes  of  computer  memory 
for  150  sound  speed  profiles. 


INTEGRATION  OF  THE  DIFFERENTIAL  EQUATIONS 

There  are  several  well-known  methods  for  integrating  a  coupled  set  of  ordinary  dif¬ 
ferential  equations  such  as  Eq.  4.  One  of  the  most  common  is  the  Runge-Kutta  method. 
Applied  to  Eq.  4,  this  method  involves  a  series  of  range  steps.  The  derivatives  of  ray 
angle,  depth,  and  travel  time  with  respect  to  range  are  calculated  at  each  step.  When 
using  any  integration  method,  it  is  important  to  monitor  truncation  error  or  to  otherwise 
assure  that  the  integration  is  accurate.  Developing  code  that  obtains  acceptable  errors  in 
the  integration  in  a  timely  fashion  is  the  essence  of  the  difficulty  of  developing  a  ray  inte¬ 
gration  code.  The  most  stringent  accuracy  requirements  occur  in  the  integration  of  ray 
angle.  The  nature  of  this  difficulty  can  be  shown  by  noting  that  the  main  term  in  Eq. 
4,(  1/c)  dc/dzcquals  d(ln(c))/dz,  and  that  sound  speed  varies  by  only  about  5%  in  the 
water  column.  Since  ln(c)  varies  by  even  less,  the  errors  in  estimating  the  derivative 
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rapidly  wreck  the  integration. 

The  Numerical  Recipes  textbook  describes  an  integration  method  involving  "adap¬ 
tive  stepping."  This  technique  involves  measuring  the  truncation  error  at  each  step  and 
then  modifying  the  step  size  such  that  the  integration  maintains  a  specified  level  of  error. 
This  ensures  both  that  the  integration  is  accurate  and  that  step  sizes  are  not  so  small  as  to 
involve  unnecessary  computation.  Properly  employed,  an  "adaptive  step"  method  can 
result  in  significant  improvements  in  computation  time. 

In  the  present  case,  however,  we  already  know  a  great  deal  about  the  ray  tracing 
problem,  unlike  a  generic  integration  problem.  For  example,  given  the  nature  of  the  ray 
tracing  problem,  it  is  probable  that  step  sizes  of  1  m  are  too  small,  and  step  sizes  of 
1000  m  are  too  large.  It  thus  appears  that  using  adaptive  stepping  will  not  result  in  an 
order  of  magnitude  decrease  in  computation  time  as  might  be  expected  for  an  unknown, 
generic  integration  problem.  It  is  true  that  adaptive  stepping  gives  a  measure  of  trunca¬ 
tion  error,  which  is  quite  important.  However,  with  experience  the  character  of  a  pre¬ 
dicted  wave  front  can  be  used  to  decide  if  predetermined  step  sizes  are  sufficiently  small, 
or  a  few  initial  calculations  can  be  made  with  different  step  sizes  to  test  the  convergence 
of  the  integration.  Once  the  parameters  of  the  integration  are  determined,  speed  of  com¬ 
putation  is  the  only  remaining  issue. 

It  is  not  presently  known  whether  an  adaptive  stepping  method  can  be  developed 
that  will  result  in  speedier  computations  than  a  user-specified  stepping.  It  is  certainly  true 
that  the  "adaptive  stepping"  method  described  in  the  Numerical  Recipes  textbook  is  more 
computationally  expensive  than  the  user-specified  stepping  method  that  has  been  devel¬ 
oped  for  the  present  ray  tracing  code.  It  may  be  that  once  the  nature  of  the  problem  has 
been  well  defined,  adaptive  step  sizing  is  a  computational  overhead  that  is  no  longer 
required.  Once  the  step  sizes  have  been  defined  for  the  problem,  the  step  sizes  do  not 
need  to  be  recalculated  every  time.  Figure  1  shows  the  step  sizes  determined  by  an  adap¬ 
tive  step  ray  trace  and  the  predetermined  step  size  currently  implemented  in  the  code. 
This  predetermined  step  size  is  linear  from  the  surface  to  a  user-specified  depth  (1500  m 
in  Figure  1)  and  has  a  hyperbolic  tangent  (tanh)  functional  form  below  that  depth.  While 
the  predetermined  steps  are  generally  considerably  smaller  than  the  adaptively  sized 
steps,  the  code  using  the  former  is  significantly  faster  than  the  code  using  the  latter. 

Besides  the  method  used  to  derive  the  step  sizes,  a  second  issue  is  the  method  of 
integration  to  be  employed.  The  Numerical  Recipes  textbook  lists  several  methods: 
Classical  4th  order  Runge-Kutta,  Cash-Karp  Runge-Kutta,  Bulirsch-Stoer,  "predictor-cor¬ 
rector,"  and  "stiff  problem"  integration.  The  Cash-Karp  Runge-Kutta  method  is  a  5th 
order  method  which  also  allows  an  estimate  of  truncation  error  using  the  difference 
between  4th  and  5th  order  Runge-Kutta  results.  Since  one  often  needs  to  integrate 
through  mixed  layers  and  other  sharp  features,  the  best  integration  method  is  probably 
Runge-Kutta,  because  many  of  the  other  methods  require  the  integration  problem  to  be 
relatively  smooth.  For  a  suitably  chosen  step  size,  the  classical  4th  order  Runge-Kutta 
and  the  Cash-Karp  Runge-Kutta  methods  produce  nearly  identical  results  in  the  ray  trace. 
As  pointed  out  in  the  Numerical  Recipes  textbook,  higher  order  does  not  necessarily 
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mean  higher  accuracy.  The  classical  4th  order  Runge-Kutta  is  much  cheaper  to  calculate. 
We  tried  to  implement  the  Bulirsch-Stoer  method,  but  it  produced  unsatisfactory  results. 
We  did  not  try  to  implement  any  "predictor-corrector"  methods. 

At  the  present  time  we  conclude  that  the  most  efficient  way  to  make  these  calcula¬ 
tions  is  to  use  the  classical  4th  order  Runge-Kutta  with  a  step  size  that  is  predetermined. 
Several  test  cases  of  varying  step  sizes  can  be  run  initially  to  ensure  that  the  integration  is 
accurate  before  proceeding  with  the  optimal  calculations.  The  main  issue  is  to  balance 
speed  and  accuracy.  We  hold  out  the  hope,  however,  that  a  more  efficient  adaptive  step¬ 
ping  scheme  may  be  devised.  (The  Bolin  RAY  code  adopts  4th  order  Runge-Kutta  and 
predetermined  step  size  as  well.  However,  the  step  sizing  is  determined  by  the  depth 
spacing  of  the  sound  speed  data  that  are  input.  As  a  result  Bowlin  RAY  predictions  using 
sound  speed  profiles  defined  by  many  depths  take  longer  to  compute  than  predictions 
using  sound  speed  profiles  defined  by  few  depths.  It  is  probably  better  to  have  the  step 
size  vs  accuracy  issue  explicitly  recognized  by  the  user  of  the  code.) 

Step  sizing  is  critical.  Unsuitably  defined  step  sizing  causes  irregular  wave  front 
predictions  and  can  introduce  biases  to  the  time  front.  The  user  performing  ray  tracing 
for  the  first  time  on  a  particular  problem  is  urged  to  try  several  step  sizing  definitions  to 
ensure  that  the  desired  accuracy  is  obtained.  It  is  far  more  difficult  to  derive  eigenrays  if 
the  time  front  is  irregular  and  rough.  For  adaptive  step  sizing  the  scaling  of  the  variables 
in  the  problem  is  critical.  High  accuracy  is  needed  at  every  step  to  avoid  the  accumula¬ 
tion  of  errors  during  the  integration.  This  means  that  the  scaling  of  the  variables  needs  to 
be  proportional  to  the  step  size.  Small  step  sizes  are  required  where  d  cos  6/ dr  ~  0  (i.e., 
where  dddz  ~  0)  or  where  the  sound  speed  changes  rapidly.  Thus  small  step  sizes  axe 
usually  required  near  the  sound  channel  axis  or  near  the  surface. 

Other  methods  of  integration  have  not  yet  been  tested.  It  may  well  be  that  faster, 
more  efficient  methods  than  Runge-Kutta  integration  could  be  found.  There  is  surely  a 
mathematical  analysis  of  the  ray  tracing  problem  that  would  point  to  the  optimal  method, 
but  this  is  beyond  the  means  of  the  authors. 


SURFACE  AND  BOTTOM  REFLECTIONS 

Surface  reflections  are  handled  in  a  simple  way  in  the  code.  The  code  determines 
at  every  step  if  a  surface  reflection  is  to  occur  in  the  next  range  step.  If  h  is  a  range  step 
to  take  after  reaching  a  depth  zn,  then  the  depth  of  the  ray  after  this  step  can  be  approxi¬ 
mated  by  zn+i  »  z„  +  dzldr  h.  If  z„+i  is  above  the  ocean  surface,  we  can  analytically 
solve  for  a  new,  shorter  step  size  such  that  the  ray  will  arrive  exactly  at  the  surface  after 
this  step.  This  works  best  (the  approximate  integration  for  the  depth  zn+\  has  less  error) 
when  very  small  step  sizes  are  taken  near  the  surface.  Once  this  has  been  done,  we  need 
merely  change  the  sign  of  the  ray  angle  and  continue  with  our  calculation. 

Because  the  step  size  derived  to  arrive  exactly  at  the  surface  is  approximate,  the 
integration  will  occasionally  produce  a  ray  that  passes  infinitesimally  above  the  ocean 
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surface.  To  ensure  that  the  calculation  can  proceed  without  mishap  in  this  case,  the  sound 
speed  table  includes  a  value  for  sound  speed  at  3  m  above  the  ocean  surface.  The  user 
can  specify  a  tolerance  for  missing  the  ocean  surface.  In  general,  depth  errors  of  several 
centimeters  or  more  in  the  surface  reflection  will  not  significantly  affect  the  ray  tracing 
results. 

Bottom  reflections  are  handled  in  a  manner  similar  to  surface  reflections,  with  a 
few  complications.  The  bottom  is  modeled  as  line  segments  between  the  depths  specified 
in  the  input  bathymetry  file,  and  the  ray  reflects  specularly.  By  simple  geometry,  if  the 
incident  angle  is  6,  the  reflected  angle  is  6 ,  and  the  angle  of  the  bottom  with  respect  to 
the  horizontal  is  <p,  then  9  =  -  0 + 2 <f>.  The  angle  <f>  is  the  arctangent  of  the  slope  of  the 
bottom. 

The  complexities  arise  in  determining  when  the  ray  is  approaching  the  ocean  bot¬ 
tom  and  in  solving  for  the  step  size  required  for  the  ray  to  reach  the  bottom  exactly  with¬ 
out  overshooting.  Suppose  the  segment  of  ocean  bottom  is  modeled  as  B(r)  =  mr  +  b, 
where  m  is  the  bottom  slope  and  b  is  a  line  intercept,  and  suppose  the  ray  path  is  z{r). 
The  sign  of  z-mr-b  therefore  determines  which  side  of  the  bottom  the  ray  is.  Similarly 
to  the  surface  reflection  procedure,  the  sign  of  ( z  +  hdz/dr)-m(r  +  h)-b  is  a  test  of 
whether  the  ray  will  cross  the  bottom  on  the  next  step  of  size  h.  This  equation  can  also 
be  used  to  solve  for  the  step  size  h  for  the  ray  to  exactly  strike  the  bottom.  However, 
because  much  larger  step  sizes  are  used  near  the  ocean  bottom  than  near  the  surface,  the 
approximate  integration  by  dzldr  does  not  work  very  well.  The  code  uses  the  constant 
ray  slope  approximation  to  determine  when  the  ray  is  near  the  bottom,  and  then  iterates 
step  size  to  obtain  the  step  size  required  to  put  the  ray  at  the  ocean  bottom  within  a  user- 
specified  tolerance.  In  practice,  a  look-up  table  of  range,  depth,  bottom  slope,  twice  bot¬ 
tom  angle,  and  the  value  of  an  intercept  of  the  line  segment  is  used  to  determine  which 
segment  of  the  bottom  the  ray  is  approaching  and  to  calculate  the  reflection  point  and  new 
ray  angle. 

This  code  does  not  model  the  finite  frequency  effects  on  ray  paths  with  small  slope 
that  pass  near  the  surface  (Murphy  and  Davis  1974;  Odom  1998).  The  geometric  ray  the¬ 
ory  described  in  this  memorandum  is  an  infinite  frequency  approximation  which  permits 
rays  to  make  an  arbitrarily  close  approach  to  the  surface  without  being  influenced  by  it. 
Only  when  the  geometric  rays  actually  contact  the  surface  are  they  reflected.  A  finite  fre¬ 
quency  ray  contains  energy  which  decays  exponentially  above  the  geometric  turning 
point,  as  shown  by  predictions  using  acoustic  modes.  The  length  of  the  exponential  tail  is 
frequency  dependent  with  lower  frequency  rays  having  longer  tails,  thus  more  energy 
above  the  geometric  turning  point.  Because  of  the  exponential  tail,  the  finite  frequency 
ray  senses  the  surface  sooner  than  the  corresponding  geometric  ray.  The  result  is  that  the 
finite  frequency  ray  turns  at  a  shorter  range  from  the  source  than  the  geometric  ray  with 
the  same  launch  angle.  This  effect  can  cause  significant  changes  to  ray  travel  times  and 
ray  paths  in  very  long  range  acoustic  transmissions.  In  one  example  from  the  Eastern 
Pacific,  the  travel  times  were  affected  by  1-2  ms  per  surface  interaction,  but  the  rays  have 
to  be  grazing  the  surface  in  just  the  right  way  for  an  effect  of  this  magnitude 
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(M.  Dzieciuch  personal  communication  1999,  Worcester  et  al.  1998).  However,  it  is  diffi¬ 
cult  to  model  the  near-surface  interaction  with  rays  in  a  general  way,  particularly  for 
broad-band  acoustic  transmissions.  For  the  time  being,  those  concerned  about  this  effect 
may  make  an  empirical  assessment  of  its  magnitude  by  comparing  ray  and  mode  or 
parabolic  equation  predictions. 


BENCHMARKING  VS  ACCURACY 

The  speed  of  a  ray  prediction  depends  on  the  accuracy  required  by  the  user.  Thus, 
it  is  virtually  impossible  to  compare  the  speeds  of  various  ray  codes,  unless  one  has  a 
means  of  ensuring  that  the  various  predictions  are  to  the  same  accuracy.  The  table  below 
shows  computation  times  for  5000  rays  when  using  an  early  version  of  the  ray  code.  The 
calculations  are  for  a  3-Mm  range,  range-dependent  sound  speed  section.  For  these  cal¬ 
culations,  the  Step  Size  Scale  is  a  scaling  factor  that  is  applied  to  the  step  size  function 
shown  in  Figure  1.  The  table  shows  that  the  computation  time  depends  linearly  on  the 
step  size. 


Step  Size  Scale  0.5  1.0  2.0  4.0 

Computation  Time  (min.)  32.0  16.5  8.75  4.5 


Figure  2  shows  that  the  time  front  predictions  associated  with  these  various  step 
size  scalings  deteriorate  considerably  as  step  size  is  increased.  By  comparison,  an  adap¬ 
tive  step  size  prediction  took  about  1  hour,  the  Bowlin  RAY  prediction  (Figure  2)  took 
34.0  minutes,  and  the  Colosi  prediction  (Figure  2)  (paranoid  for  accuracy)  took  over 
5  hours.  While  for  all  of  these  predictions  the  location  of  the  individual  ray  arrivals  may 
vary,  all  of  the  time  fronts  agree  within  milliseconds.  The  shape  and  absolute  travel  times 
of  time  fronts  appear  to  be  fairly  robust. 

The  lesson  here  is  that  the  user  of  a  ray  code  needs  to  balance  his  or  her  desire  for 
accuracy  against  the  speed  of  computation. 

The  astute  reader  will  have  noticed  in  Figure  2  occasional  gaps  in  the  time  fronts, 
particularly  in  the  earliest  part  of  the  arrival  pattern.  These  gaps  apparently  result  from 
the  ray  approximation;  similar  predictions  using  the  parabolic  equation  do  not  show  these 
gaps.  In  the  ray  approximation,  rays  which  pass  near  a  surface  layer  will  either  travel  into 
that  surface  layer  or  miss  it  altogether,  even  though  the  launch  angle  in  the  two  cases  may 
differ  only  infinitesimally.  This  property  of  the  ray  approximation  is  a  likely  cause  of  the 
gaps  in  the  time  front  prediction.  When  such  gaps  occur  near  the  depth  of  a  receiver,  it 
can  be  difficult  to  obtain  an  associated  eigenray  prediction.  Indeed,  difficulty  in  obtaining 
the  eigenray  prediction  for  a  ray  known  to  exist  probably  results  from  such  an  unphysical 
gap  in  the  time  front. 
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CONCLUSIONS 

We  have  described  a  method  for  making  fast,  accurate  computations  of  acoustic 
rays  as  a  tool  for  ocean  acoustic  tomography.  We  described  ray  equations,  including  a 
derivation  of  ray  equations  that  include  the  effects  of  weak  ocean  currents  (Appendix  F). 
Methods  using  cubic  spline  interpolation  and  a  look-up  table  allow  sound  speed  and 
sound  speed  gradient  to  be  calculated  rapidly  and  accurately  at  arbitrary  range  and  depth. 
The  choice  of  the  step  size  used  in  integrating  the  differential  equations  is  critical,  affect¬ 
ing  both  the  computation  time  and  the  accuracy  of  the  ray  predictions.  The  best  method 
for  integrating  the  ray  equations  appears  to  be  ordinary,  classical  4th  order  Runge-Kutta 
integration.  Presently  known  integration  methods  using  the  adaptive-step  techniques  to 
maintain  a  user-specified  numerical  accuracy  appear  to  carry  too  much  computational 
cost  to  be  competitive.  A  user-specified  step  size  greatly  increases  the  efficiency  of  the 
computation,  but  a  more  efficient  adaptive-stepping  method  may  yet  be  devised.  Integra¬ 
tion  accuracy  can  be  initially  checked  for  a  particular  problem  by  repeating  the  ray  pre¬ 
dictions  with  several  trial  step-size  functions. 

For  5-Mm  range  acoustic  transmissions  across  the  North  Pacific,  eigenray  travel 
times  calculated  with  the  code  described  here  agreed  with  the  travel  times  of  the  Bowlin 
code  to  within  20  ms.  Similar  agreement  was  found  for  other  calculations  at  similar  mul¬ 
timegameter  range.  The  upper  and  lower  turning  depths  of  eigenray  paths  calculated 
using  the  code  here  and  the  Bowlin  code  agreed  to  within  a  few  meters.  Sound  speeds 
calculated  using  the  annual  mean  Levitus  ’94  ocean  atlas  were  used  for  these  calcula¬ 
tions.  When  calculating  time  fronts  to  similar  accuracy,  the  present  code  appears  to  be 
comparable  in  speed  to  the  Bowlin  code.  However,  when  calculating  eigenrays 
(Appendix  B),  the  present  code  appears  to  be  2-3  times  faster  than  the  Bowlin  code 
because  most  of  the  time  it  integrates  only  two  differential  equations.  In  addition,  the 
present  eigenray  code  appears  to  be  a  little  more  effective  at  finding  the  eigenrays  than 
the  Bowlin  code,  so  fewer  rays  need  to  be  calculated  to  define  the  initial  time  front  used 
to  find  the  eigenrays. 

The  code  described  here  is  inherently  unstable,  yet  it  is  highly  flexible.  The  com¬ 
puter  code  is  easily  modified  for  particular  problems.  One  application  of  this  code  may 
be  to  calculate  the  forward  problem  matrices  used  for  travel  time  inversion  while  calculat¬ 
ing  the  eigenrays.  This  calculation  would  result  in  more  accurate  matrices  and  it  will  not 
require  saving  the  ray  paths  to  one’s  hard  disk.  As  described  in  some  of  the  appendices, 
the  code  can  calculate  eigenrays,  can  be  modified  to  run  on  a  parallel  computer  for 
rapidly  obtaining  ray  predictions,  or  can  be  modified  to  include  the  effects  of  ocean  cur¬ 
rent.  The  present  suite  of  software  consists  of  about  1800  lines  (including  lots  of  com¬ 
mentary)  of  FORTRAN  code. 
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APPENDIX  A 

A  TECHNICAL  SUMMARY  AND  A  FLOW  CHART 
OF  THE  COMPUTER  CODE 

A  flow  chart  sketching  the  computer  code  is  shown  in  Figure  3.  The  main  program 
ray.f  first  sets  up  the  sound  speed  look-up  tables  (subroutine  speed.f)  before  integrating 
the  rays  (subroutine  dodeint.f)  with  the  desired  ray  launch  angles.  The  code  is  in  double 
precision  throughout.  Note  that  for  the  code  to  be  universally  portable  the  common  tables 
must  have  the  double  precision  arrays  and  variables  listed  before  the  integer  variables. 

The  sound  speed  look-up  tables  (ctab)  are  calculated  using  the  Numerical  Recipes 
cubic  spline  routines  dspline.f  and  dsplintf.  The  subroutine  dspline.f  calculates  the  sec¬ 
ond  derivative  of  sound  speed,  from  which  the  cubic  splined  sound  speeds  are  obtained. 
The  subroutine  dsplint.f  has  been  modified  to  obtain  the  subroutine  dsplint_both.f  which 
returns  the  values  of  both  sound  speed  and  sound  speed  gradient.  While  reading  in  the 
values  for  sound  speed,  the  subroutine  also  obtains  the  table  of  profile  ranges  (range). 

A  step  size  look-up  table  (steps)  is  also  determined  in  the  subroutine  speed.f.  As 
mentioned  earlier,  this  predetermined  step  size  is  linear  from  the  surface  to  a  user-speci¬ 
fied  depth  and  then  has  a  tanh  functional  form  below  that  depth.  The  user  is  asked  to 
specify  the  transition  depth  from  linear  to  tanh  forms  and  the  step  size  values  at  the  sur¬ 
face,  at  the  transition  depth,  and  at  5500-m  depth. 

The  bathymetry  look-up  table  (btab)  is  determined  in  the  subroutine  bathy.f.  This 
subroutine  reads  in  the  maximum  allowable  bottom  bounces  (MBONK),  the  tolerance  for 
the  ray  to  miss  the  ocean  bottom  (btol),  and  the  bathymetry  data.  The  bathymetry  look¬ 
up  table  consists  of  the  five  variables:  range,  depth,  slope,  twice  the  slope  angle,  and  a 
value  of  an  intercept  of  the  line  segment. 

Once  the  sound  speed,  range,  step  size  and  bathymetry  look-up  tables  have  been 
defined,  the  rays  are  integrated  using  the  subroutine  dodeint.f.  Presently  the  code  loops 
over  a  set  of  rays  at  equally  incremented  launch  angles.  The  subroutine  dodeint.f  has 
been  heavily  modified  from  the  Numerical  Recipes  subroutine  example  of  the  same  name. 
The  role  of  this  subroutine  is  as  a  driver  for  taking  the  Runge-Kutta  steps  (step  size  h) 
with  the  subroutine  drk4.f;  the  user  wishing  to  change  the  method  of  integration  from 
Runge-Kutta  may  merely  substitute  an  alternate  subroutine  to  drk4.f.  The  subroutine 
drk4.f  has  been  only  slightly  modified  from  the  Numerical  Recipes  subroutine  of  the 
same  name.  The  subroutine  drk4.f  calls  the  subroutine  derivs.f  which  calculates  the  cou¬ 
pled  differential  equations  (the  ray  equations)  at  arbitrary  depth  and  ray  angle.  The  sub¬ 
routine  dodeint.f  saves  the  ray  paths  derived  during  the  integration  in  variables  xp  and  yp 
and  returns  the  cosine  and  sine  ray  angle  at  the  receiver  range,  the  depth  of  the  ray  at  the 
receiver  range,  and  the  ray  travel  time  for  writing  to  a  file. 

Range  dependence  is  implemented  in  the  derivs.f  subroutine  using  the  variable 
deltaX;  deltaX  is  the  horizontal  range  between  a  step  in  the  integration  and  the  range  of 
the  sound  speed  profile  most  recently  exceeded  by  the  ray.  Since  range  dependence  is 
implemented  using  a  constant  sound  speed  gradient  between  profiles,  the  sound  speed  at 
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any  step  in  the  integration  can  be  found  by  c(r,  z )  =  c(rh  z)  +  deltaX  dc/dr.  Here, 
c(r,,  z)  is  the  sound  speed  profile  at  range  r,  .  The  value  for  the  vertical  gradient  of  sound 
speed  can  be  found  similarly.  In  general,  the  step  size  h,  0(5-200)  m,  is  much  less  than 
deltaX,  0(10  km). 
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APPENDIX  B 

CALCULATION  OF  EIGENRAYS 

There  are  several  ways  in  which  eigenrays  can  be  computed  accurately,  yet  with  far 
fewer  calculations  than  are  required  for  a  complete  time  front  prediction.  In  all  cases, 
eigenrays  must  be  found  by  first  tracing  a  fan  of  rays  and  then  determining  the  ray  launch 
angles  that  result  in  rays  arriving  at  the  receiver  depth.  However,  there  are  a  number  of 
shortcuts  that  can  be  implemented. 

In  a  ray  fan  equally  spaced  in  launch  angle,  most  of  the  rays  concentrate  in  the 
cusps  of  the  time  front  and  so  frequently  do  little  to  resolve  the  time  front  in  the  depth 
region  of  interest.  Thus,  one  of  the  simplest  ways  to  speed  up  eigenray  calculations  is  to 
use  an  initial  prediction  with  a  small  number  of  rays  in  order  to  define  the  range  of  ray 
angles  that  arrive  near  the  depth  of  the  receiver.  Using  these  selected  angle  ranges,  a  sec¬ 
ond  prediction  with  many  rays  can  be  done  that  is  far  more  efficient  at  defining  the  time 
front  near  the  depth  of  the  reciever. 

A  second  way  to  gain  efficiency  is  to  omit  the  calculation  of  travel  time  in  the  ini¬ 
tial  fan  of  rays.  All  that  is  required  to  derive  eigenrays  is  the  angle  that  produces  a  ray 
that  arrives  at  the  receiver  depth,  so  that  travel  time  is  not  necessary.  The  present  incarna¬ 
tion  of  the  eigenray  code  implements  this  using  only  two  coupled  differential  equations 
(for  tangent  of  ray  angle  and  ray  depth;  see  Eq.  B1  below),  rather  than  the  four  required 
for  obtaining  travel  time  as  well.  Once  the  launch  angles  resulting  in  rays  that  arrive  at 
the  receiver  depth  are  found,  a  small  number  of  integrations  with  all  four  differential 
equations  can  then  be  performed  to  obtain  all  the  information  about  the  eigenrays. 

Eigenray  predictions  are  probably  more  efficient  if  fewer,  more  accurate  rays  are 
used  to  calculate  the  time  front  rather  than  many,  inaccurate  rays.  The  more  linear  the 
relation  between  launch  angle  and  arrival  depth,  the  easier  it  will  be  to  obtain  the  eigen¬ 
rays. 

After  the  two  ray  arrivals  that  span  the  receiver  depth  are  determined,  eigenrays  are 
probably  best  found  by  a  sequence  of  bisections  of  ray  launch  angle  giving  a  sequence  of 
rays  that  converge  on  the  receiver  depth.  An  initial  interpolation  to  solve  for  the  eigenray 
launch  angle  might  improve  efficiency,  but  in  general  the  ray  arrivals  are  not  linear  in 
launch  angle  and  arrival  depth,  so  the  interpolation  will  not  be  accurate.  The  code 
presently  uses  ordinary  bisection  on  rays  that  span  the  receiver  depth.  Because  of  non- 
linearities  in  the  relation  between  launch  angle  and  receiver  depth,  the  sequence  of  bisec¬ 
tions  sometimes  misses  the  receiver  depth.  Various  contingencies  have  been  built  into  the 
code  to  account  for  this,  but  it  may  sometimes  happen  that  the  eigenray  cannot  be  found 
at  all.  It  appears  that  at  3-Mm  range,  usually  less  than  five  bisections  are  required  to  find 
a  ray  that  arrives  within  10  m  of  the  receiver  depth. 

If  one  wants  only  the  travel  time  of  the  eigenrays,  one  could  interpolate  the  pre¬ 
dicted  time  front  to  get  the  travel  times  associated  with  the  depth  of  the  receiver.  Travel 
times  determined  in  this  way  will  (for  benign  oceans)  be  accurate  to  within  a  few  mil¬ 
liseconds. 
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Thus,  the  travel  times  of  the  eigenrays  can  be  determined  without  having  to  actually  solve 
for  the  eigenray  path.  This  technique  has  not  been  implemented. 

The  code  here  is  flexible  enough  that  the  user  can  easily  implement  ideas  for  find¬ 
ing  eigenrays.  If  someone  develops  code  that  works  well  for  determining  eigenrays,  he  or 
she  is  requested  to  forward  those  methods  to  the  authors  of  this  report. 

The  code  for  finding  eigenrays  is  similar  to  that  described  in  Appendix  A,  with  the 
addition  of  subroutines  dodeint_short.f,  drk4_short.f,  and  derivs_short.f.  These  sub¬ 
routines  use  only  two  differential  equations  in  tangent  of  ray  angle  and  ray  depth.  If 
tn  =  tan(6),  then  these  equations  are 


dtfl  /,  2x  dz-c 

—  =  (l  +  m2)-^ 
dr  c 


(Bla) 


dz 


(Bib) 


dt  _  Vl  +  m2^ 
dr  c 

J 


(Blc) 


The  subroutine  used  to  find  the  eigenray  ray  launch  angles  is  find_arr.f.  This  rou¬ 
tine  uses  the  time  front  information  from  the  initial  ray  traces  to  select  pairs  of  ray  angles 
with  depths  at  the  receiver  range  that  span  the  receiver  depth.  A  series  of  bisections  on 
ray  launch  angle  is  used  to  converge  on  the  desired  eigenray  ray  launch  angle.  The  rays 
at  the  selected  angles  are  then  recalculated  using  all  four  differential  equations  to  obtain 
all  the  information  about  the  ray.  Since  this  entire  procedure  relies  on  these  two  different 
integrations  giving  nearly  identical  ray  paths,  the  procedure  is  actually  one  test  of  integra¬ 
tion  accuracy.  If  the  eigenrays  in  the  second  integration  miss  the  receiver  depth,  one 
problem  may  be  that  the  step  sizes  are  too  big.  However,  the  integrations  are  so  sensitive 
that  the  series  of  step  sizes  during  the  two  integrations  must  be  fairly  well  synchronized. 
If  even  slightly  different  step  sizes  are  used  in  the  second  integration,  the  rays  can  miss 
the  the  receiver  depth  by  0(  100  m). 

Once  the  eigenrays  have  been  calculated,  frequently  the  next  step  for  the  tomogra- 
pher  is  to  calculate  the  forward  problem  matrix, 


_  r  Tjj(x)ds 
iJ  J  c(x)2 

1  i 


(B2) 
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Here  T,  is  the  ith  ray  path,  t}j  is  one  of  a  set  of  functions  used  to  model  the  ocean,  ds  is 
an  element  of  ray  path  length  (equal  to  dr/cos(0)),  and  c  is  sound  speed.  Traditionally,  to 
facilitate  this  calculation  the  parameter  ds/c2  is  calculated  for  each  saved  step  in  the  ray 
trace.  With  this  parameter,  the  above  integral  can  be  calculated  with  relative  ease.  How¬ 
ever,  ds/c2  is  not  entirely  easy  to  calculate  with  rigorous  accuracy. 

The  problem  of  calculating  ds/c2  is  illustrated  in  Figure  4.  While  a  ray  path  is  cal¬ 
culated  using  very  small  horizontal  step  sizes,  0(5—200  m),  the  ray  path  is  saved  at  a 
much  large  range  increment,  0(1000  m).  What  is  required,  then,  is  a  value  for  ds/c 2 
appropriate  for  this  much  larger  horizontal  range  increment.  It  is  not  obvious  how  this 
parameter  is  to  be  calculated  rigorously.  It  is  easy,  of  course,  to  approximate  ds  by  sum¬ 
ming  the  small  dsj  that  occur  during  the  integration,  but  what  of  1/c2  ?  Should  the  values 
of  c,  1/c,  or  1/c2  be  averaged  over  the  several  small  steps?  Sound  speed  varies  so  little 
over  a  path  increment  that  it  probably  does  not  matter  too  much.  The  present  code  calcu¬ 
lates  dsj/c2  at  each  integration  step  and  sums  these  values  over  the  several  integration 
steps  that  comprise  the  saved  range  increment.  The  value  of  ds/c 2  is  saved  together  with 
the  ray  path  range  and  depth  pairs.  Possible  errors  may  arise  at  depths  where  the  sound 
speed  or  mode  functions  are  varying  rapidly,  however. 

With  the  present  set  of  subroutines,  it  would  be  relatively  easy  to  merge  the  ray 
tracing  routines  with  the  calculation  of  the  forward  problem.  In  this  way,  the  forward 
problem  matrices  can  be  calculated  at  the  very  small  step  sizes  of  the  integration  on  the 
fly,  as  it  were,  without  having  to  fill  one’s  hard  disk  up  with  the  stored  ray  paths.  This  is 
the  ideal  way  to  calculate  the  forward  problem  matrices. 

The  Bowlin  et  al.  (1992)  RAY  code  appears  to  address  the  problem  of  ds/c2  by 
integrating  a  variable  q  =  s/c 2 .  With  this  variable,  ds/c 2  may  presumably  be  calculated  at 
the  nth  step  by  qn  -  qn-\.  This  procedure  is  not  correct,  since  d(s/c2)  =  ds/c2  -2s dc/c 3. 
The  second  term  has  a  dangerous  dependence  on  the  path  length  itself.  The  second  term 
oscillates  in  sign  because  dc  changes  sign  from  upward-going  to  downward-going  rays. 
Calculation  of  this  variable  requires  the  integration  of  an  additional  differential  equation. 
For  small  values  of  s,  the  values  of  ds/c2  derived  from  present  code  and  from  the  Bowlin 
code  appear  agree  to  1  part  in  10,000  for  the  North  Pacific  environment  modeled  by  the 
annual  mean  Levitus  atlas. 
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APPENDIX  C 

MODIFYING  THE  CODE  FOR  COMPUTATIONS  IN  PARALLEL 

Calculating  ray  predictions  for  long-range  transmissions  is  a  natural  application  for 
parallel  computation.  Each  ray  trace  is  independent  of  all  the  others,  so  several  computer 
processors  can  calculate  subsets  of  rays  separately.  A  parallel  version  of  the  ray  code 
described  here  has  been  developed  using  the  MPI  suite  of  subroutines  (Gropp  et  al. 
1994). 

In  order  for  several  processors  to  be  able  to  trace  rays,  all  that  is  required  is  that 
they  have  access  to  the  sound  speed  look-up  table.  Thus,  the  parallel  version  of  the  code 
broadcasts  this  look-up  table,  together  with  a  few  other  variables,  to  all  available  proces¬ 
sors.  With  this  information,  the  available  processors  can  trace  subsets  of  the  desired 
rays.  When  they  have  completed  these  calculations,  they  then  communicate  the  results 
(travel  times,  ray  paths,  etc.)  back  to  the  master  process  for  printing  or  other  use.  The 
MPI  routines  perform  these  communications  tasks  (broadcasting  to  all  processors  or 
sending  and  receiving  information  between  processors)  relatively  easily.  The  speed  of 
the  ray  calculations  scales  well  with  the  number  of  available  processors. 

With  a  homogeneous  set  of  processors,  the  complete  set  of  rays  to  be  calculated 
can  be  broken  up  into  equally  sized  subsets.  Each  processor  then  gets  a  subset  of  rays  to 
calculate.  The  present  parallelized  version  of  the  ray  code  adopts  this  strategy. 

With  an  inhomogeneous  set  of  processors,  "load  balancing"  becomes  an  issue. 
Fewer  rays  should  be  sent  to  slower  processors.  It  may  be  better  to  assign  rays  singly  to 
processors  and  assign  additional  rays  to  a  particular  processor  as  its  ray  calculations  are 
completed.  MPI  includes  "non-blocking"  send  and  receive  communication  subroutines 
which  may  be  used  for  this  purpose.  That  is,  the  code  will  not  wait  for  a  particular  ray 
trace  to  be  completed  before  assigning  rays  to  other  processors.  In  this  way,  a  processor 
that  is  four  times  faster  than  another  processor  will  complete  four  times  as  many  rays 
without  requiring  the  user  to  preallocate  unequally  sized  subsets  of  rays  to  the  processors. 
This  strategy  has  not  been  implemented,  but  it  would  not  be  difficult  to  do  so. 
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APPENDIX  D 

INPUT  AND  OUTPUT  FILES  AND  OTHER  OPERATIONAL 

INFORMATION 


As  noted  often  earlier,  this  code  is  meant  to  be  modified  by  the  user  for  his  or  her 
own  purposes.  This  includes  the  formats  of  the  sound  speed  data  files  and  the  means  to 
input  other  information  such  as  receiver  depth,  range,  and  other  parameters  of  the  ray 
trace  as  well  as  the  output  formats.  The  present  setup  of  the  code  is  described  here  as  an 
example  of  a  way  to  implement  the  code. 

The  code  is  presently  run  by  executing  the  command 

time  eigenray  <  in.ray 

"time",  of  course,  reports  the  time  it  takes  to  execute  the  raytrace;  "eigenray"  is  the 
executable,  and  the  file  "in.ray"  contains  all  of  the  relevant  parameters  that  "eigenray" 
asks  for. 


An  example  of  a  file  for  input  to  eigenray: 


8.  12. 

100  1500. 


100 

1000  1000 
0.0  3000000 
1000.0 
10.0 
1 
0 
1 

test.ssp 

1500. 

5.  100.  500. 

1.0 

30 

5.0 

test.bth 


Starting  and  ending  values  for  ray  angle  (degrees) 

#  rays  to  sketch  out  time  front,  and  depth  below  which  to  omit  rays  in  the 
second  raytrace.  The  second  ray  trace  is  not  performed  for  time  front 
calculations. 

#  rays  in  second  ray  trace. 

Source  and  receiver  depths  (meters) 

Start  and  end  ranges  (meters) 

The  range  increment  to  save  the  ray  paths  (meters) 

Tolerance  for  missing  receiver  depth  (meters) 

Save  eigenray  paths?  0=no,  l=yes 

Calculate  either  eigenrays  (=0)  or  a  time  front  (=1) 

Enable  ocean  bottom?  0=no,  l=yes 
The  sound  speed  filename 

Transistion  depth  for  step  size,  linear  to  tanh  (meters) 

Values  of  step  size  at  the  surface,  the  transition  depth,  and  at  the  bottom 
(all  in  meters) 

Step  size  scaling  e.g.,  0.5, 1.0, 2.0;  for  testing  accuracy 

The  maximum  number  of  bottom  bounces  before  dropping  a  ray 

The  tolerance  for  missing  the  ocean  bottom  (meters) 

The  bathymetry  file  name.  These  last  three  lines  are  not  required 
if  the  flag  disabling  bathymetry  is  set. 
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The  file  of  sound  speeds  has  the  same  format  as  that  of  the  Bowlin  RAY  code,  with  "-1 
range”  marking  the  start  of  a  profile,  followed  by  depth  and  sound  speed. 

An  example  sound  speed  file: 

-1  0.000 
0  1500.597 
10  1500.465 
20  1499.793 
30  1498.795 
50  1495.540 
75  1491.698 
100  1489.109 
125  1487.366 
150  1486.174 
200  1484.450 
250  1483.021 
300  1481.752 
400  1479.939 
500  1479.103 
600  1479.168 
700  1479.679 
800  1479.737 
900  1480.225 
1000  1480.835 
1100  1481.586 
1200  1482.276 
1300  1483.210 
1400  1484.137 
1500  1485.137 
1750  1487.678 
2000  1491.187 
2500  1498.126 
3000  1506.161 
3500  1514.534 
4000  1523.252 
4500  1532.259 
5000  1541.404 
5500  1550.529 

-1  50.963  A  new  sound  speed  profile  at  50.963  km  range. 

0  1500.751 
10  1500.612 
20  1499.982 

30  1499.038  (etc.,  repeated  as  necessary) 
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When  the  flag  enabling  bathymetry  is  set,  the  code  requires  a  bathymetry  file  consisting 
of  simply  two  columns  of  range  and  depth,  both  in  meters. 

The  time  front  or  eigenray  data  are  written  to  the  file  "ray.info".  This  file  consists 
of  the  travel  time,  ray  angle  at  the  source,  ray  angle  at  the  receiver,  ray  upper  and  lower 
turning  depths,  depth  of  ray  at  receiver  range,  and  number  of  turning  points.  The  ray 
paths  are  written  to  the  file  "ray.data";  this  is  not  recommended  when  calculating  a  time 
front.  The  columns  of  this  file  are  range  (meters)  and  depth  (meters)  along  the  ray  path 
and  the  parameter  ds/c 2  along  the  ray  path  which  is  used  for  calculating  the  matrices  for 
inversion  of  travel  times. 
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APPENDIX  E 

THE  RAY  EQUATIONS  IN  TERMS  OF  SOUND  SLOWNESS 

The  ray  equations  (4a-4d)  can  be  considerably  simplified  if  they  are  formulated  in 
terms  of  sound  slowness.  Sound  slowness  is  the  reciprocal  of  sound  speed,  p  =  1/c.  The 
following  equations  are  equivalent  to  Eqns.  (4a-4d),  but  they  are  not  currently  imple¬ 
mented  in  the  code.  They  are  included  here  for  reference.  If  A  =  p  cos{8 )  and 
B  —  p  sin{6),  then 


dA  _  p  dp 
dr  A  dr 


(Ela) 


dB  p  dp 
dr  A  dz 


(Elb) 


dz_B 
dr  A 


(Elc) 


dt  p2 
dr  A 


(Eld) 


These  equations  are  ordinary  differential  equations  for  A,  B,  z,  and  t.  One  advantage  of 
using  these  equations  is  that  one  division  is  saved,  and  so  these  equations  may  be  about 
15%  faster  to  compute  than  Eqns.  (4a-4d).  However,  in  Eqns.  (4a— 4d),  the  terms  involv¬ 
ing  dc/dr  are  neglected,  while  this  is  a  term  essential  to  equation  (Ela).  It  is  unclear  to  us 
what  the  implications  of  this  fact  are,  e.g.,  whether  small  errors  in  the  calculation  of  dc/dr 
may  prove  to  be  a  difficulty  in  implementing  Eqns.  (Ela-Eld).  The  approximation  that 
dA/dr  =  0  is  probably  not  adequate.  In  addition,  there  are  certain  logistical  difficulties 
with  using  a  look-up  table  for  sound  slowness. 


TM  3-98  21 


UNIVERSITY  OF  WASHINGTON  ‘APPLIED  PHYSICS  LABORATORY 


APPENDIX  F 

THE  RAY  EQUATIONS  WITH  CURRENT 

The  code  does  not  presently  offer  the  option  of  including  the  effects  of  ocean  cur¬ 
rent.  We  derive  here  the  ray  equations  that  include  the  effects  of  ocean  current.  The  user 
interested  in  calculating  the  effects  of  current  can  modify  the  code  to  include  the  current 
effects. 

A  quick  review  of  the  literature  (i.e.,  Franchi  and  Jacobson  1972;  Munk  et  al.  1995) 
found  no  ray  equations  that  included  the  effects  of  current  and  could  be  easily  imple¬ 
mented  for  numerical  calculation.  There  seems  to  be  a  tendency  in  the  literature  toward 
the  cryptic  and  toward  obscure  notation.  We  therefore  include  here  a  derivation  which 
may  prove  to  be  practically  useful. 

The  textbook  by  Pierce  (1989)  gives  ray  equations  with  current  in  terms  of  the 
derivative  with  respect  to  travel  time.  These  equations  are 


dr 

~dt 


(Fla) 


—  =  -  —  Vc  -  p  x  (V  x  v)  -  (p  •  V)v 
dt  c 


or  (in  Cartesian  coordinates) 


(Fib) 


dpi  _  Q  5c  _  3  a 
dt  C  drt  hPj  dr,  j 


(F2) 


Here  r  is  ray  spatial  coordinate,  p  =  k lco  is  wave-slowness  vector,  v  is  current,  and 

£2  = - —  (F3) 

c  +  v  •  n 

where  n  is  a  unit  normal  that  is  parallel  to  the  ray.  We  can  immediately  make  a  number 
of  simplifications  to  these  equations.  We  assume  that  v  =  (v,0,0),  i.e.,  current  is  only  in 
the  direction  of  the  plane  of  acoustic  propagation.  We  assume  that  the  user  of  the  code 
will  project  current  onto  the  path  between  the  source  and  receiver  and  that  currents  in  the 
vertical  are  negligible.  Further,  the  wave-slowness  vector  is 


p  =  {cos  6,  0,  sin  6)1  c 


(F4) 
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and 


Q  = 


C  +  V  cos  0 

These  simplifications  result  in  the  following  equations: 

dr  c  cos  6 


dt  Q 

dz  _c  sin  0 
dt  Q 


+  v 


(F5) 


(F6a) 


(F6b) 


d 


(  cos  0  \ 


d 


k  c  J 

3c 

cos0  3v 

dt 

c  3r 

c  dr 

(  sin  0  \ 

l  C  J 

£2  3c 

cos0  3v 

dt 

c  3  z 

c  3  z 

(F6c) 


(F6d) 


While  it  may  make  for  an  interesting  exercise  to  integrate  these  equations  using 
small  steps  of  travel  time,  this  is  not  particularly  useful  in  the  present  context.  We  there¬ 
fore  seek  to  use  the  magic  of  mathematics  to  transform  these  equations  into  something 
useful.  We  will  assume  that  v/c  «  1  and  keep  terms  only  to  first  order  in  v/c.  For  exam¬ 
ple, 


Q  =  1  - 


VCOS0 

c 


(F7) 


The  equation  for  dtldr  can  be  found  by  taking  the  reciprocal  of  the  equation  for  dr/dt: 


dt 

dr 


1 


c  cos  0 


(2  +  tan2  0) 


(F8a) 


The  equation  for  dzldr  can  be  similarly  found  by  dzldtldrldt : 


TM  3-98  23 


UNIVERSITY  OF  WASHINGTON  •  AP PLI ED  PH YSICS  LABORATORY 


dz 

dr 


—  tan  6\  1 


V  A 

c  cos  6  ) 


(F8b) 


The  differential  equations  for  d  cos  61 dr  and  d  sin  61  dr  can  be  derived  similarly  by  noting 
that 


( cos  6  \ 


v  £  ; 

dr 


1  dcos6 
c  dr 


1 

—  cos  6 
c 2 


dc 

dr 


(F9) 


and 


dc  dc  dz  dc 
dr  dr  dr  dz 


(F10) 


After  six  pages  of  algebra  (which  interested  readers  should  probably  do  for  themselves  to 
check),  the  set  of  differential  equations,  including  those  for  ray  angles,  is 


d  cos  6 

Jr 


=  -  sin  6\ 


( 1  dc 


1  dc 


—  tan  6  - 
c  dr  c  dz 


(FI  la) 


+  (3  +  tan2  6) 


v  dc 
c2  dr 


1  dv 
c  dr 


v  dc 

tmesTz 


d  sin  6 

Jr 


■  cos  6\ 


1  dc  1  dc} 

cTrt:me--cTz  ) 


^  v  dc 


1  dv 
c  dz 


(FI  lb) 


dz 

dr 


=  tan  6  1  - 


— 1 
CCOS0  J 


(Flic) 


—  =  — - 4  (2  +  tan2  6)  (Flld) 

dr  c  cos  6  c2 

Equations  El  la-El  Id  constitute  the  ray  equations,  modified  to  account  for  ocean  current 
to  first  order.  Note  that  when  v  =  0  these  equations  satisfactorily  reduce  to  the  original 
ray  equations,  Eqs.  4a-4d. 
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In  practice,  it  would  be  easy  to  include  the  current  terms  in  the  subroutine  derivs.f 
defining  the  differential  equations.  The  second  modification  would  be  in  the  subroutine 
speed.f  to  create  a  look-up  table  for  current  analogous  to  the  table  for  sound  speed. 
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Step  Size  (m)  Step  Size  (m) 

Figure  1 .  Step  sizes  determined  by  an  adaptive  ray  trace  and  the 
predetermined  step  size  presently  implemented  in  the  code.  The  left 
panel  shows  the  step  sizes  as  a  function  of  depthderived  by  an  adaptive 
step  size  algorithm  used  to  calculate  a  single  ray  with  a  fairly  small  error 
tolerance.  The  right  panel  shows  an  analytic  step  size  function  that 
produces  ray  predictions  as  accurate  as  the  adaptive  step  sizing  but  that 
require  much  less  computation  time. 
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Figure  2.  Time  front  predictions 
associated  with  various  step  size 
scalings.  The  left  panels  show  the 
deterioration  of  the  time  front  as  the 
step  size  is  increased.  For  very  small 
step  sizes  the  ray  arrivals  are  evenly 
spaced.  For  comparison,  predictions 
using  the  Bowlin  and  Colosi  codes  are 
shown  in  the  right  panels.  In  all  cases 
the  shape  and  travel  time  of  the  time 
front  are  nearly  identical. 
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Program  EIGENRAY 


Information  about 
source,  receiver 

depth,  - ► 

range,  etc. 


Subroutine  SPEED 

step  size  table  -  steps(#depths) 

(read  in  steps  parameters) 

sound  speed  table  - 

ctab(#depths,#profiles,6) 

(read  in  sound  speed  profiles) 


Subroutine  Subroutine 
SPLINE  SPLINT_BOTH 
(2nd  Deriv.)  (c(z),  dc/dz) 


Loop  over 
ray  launch 
angles 


* 


> 


Write  out  ray 
information 


Stop 


Subroutine  ODEINT 
integrate  over  a  path 
(save  information  at  regular 
intervals  if  desired) 


One  call  to  A  multitude 

initialize  of  steps 


Subroutine  Subroutine 

DERIVS  RK4 

Runge-Kutta 

step 


Figure  3.  A  flow  chart  of  one  implementation  of  the 
ray  trace  code.  This  code  is  designed  to  be  easily 
modified  for  the  user?s  purposes. 


Subroutine 
DERIVS 
(ray  eqns.) 
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Depth 


c(z)  or  ri(z) 


Saved  range  increment  =  E  hj 


Range 


Figure  4.  A  schematic  figure  showing  the  relation  of  increments 
of  ray  path,  s,  to  the  sound  speed  profile  or  mode  function.  The 
ray  path  is  calculated  at  small  range  increments  hj  (e.g.,  5-200  m) 
but  saved  at  a  much  larger  increments  (e.g.,  1000  m).  A  value 
of  ds/c2  associated  with  the  saved  path  increment  is  required 
to  calculate  the  matrices  used  in  the  inversion  of  travel  times. 

It  is  probably  better  to  calculate  the  matrices  while  performing  a 
ray  trace;  better  accuracy  may  be  obtained  by  using  the  small 
step  sizes. 
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